package com.jellyfishumbrella.noise;

import java.util.Random;

import com.jellyfishumbrella.GL.Vector3d;

public class PerlinNoise {
    private static short permute[];
    private static short grad[];
    private static Vector3d vect[];

    private static final int TABLE_SIZE = 256;
    private static final double GRAD_SCALE = 2.0 / Math.sqrt(2.0);

    static {
        int i, j;
        short temp;
        Random random = new Random();
        Vector3d v;
        double len;

        // Precalculate the permutation table.

        random.setSeed(0);
        permute = new short[TABLE_SIZE * 2];
        grad = new short[TABLE_SIZE * 2];
        for (i = 0; i < TABLE_SIZE; i++)
            permute[i] = (short) i;
        for (i = 0; i < TABLE_SIZE; i++) {
            j = Math.abs(random.nextInt()) % TABLE_SIZE;
            temp = permute[i];
            permute[i] = permute[j];
            permute[j] = temp;
        }
        for (i = 0; i < TABLE_SIZE; i++)
            grad[i] = (short) (permute[i] % 12);

        // Precalculate the random vectors for the vector noise function.

        vect = new Vector3d[TABLE_SIZE * 2];
        for (i = 0; i < TABLE_SIZE; i++) {
            do {
                v = new Vector3d(2.0 * random.nextDouble() - 1.0,
                        2.0 * random.nextDouble() - 1.0, 2.0 * random
                                .nextDouble() - 1.0);
                len = v.distance(v);
            } while (len > 1.0);
            v.scale(1.0 / Math.sqrt(len));
            vect[i] = v;
        }

        // Extend all the tables to be twice as long.

        for (i = 0; i < TABLE_SIZE; i++) {
            permute[TABLE_SIZE + i] = permute[i];
            grad[TABLE_SIZE + i] = grad[i];
            vect[TABLE_SIZE + i] = vect[i];
        }
    }

    /** Given a point in 3D space, return the value of the scalar noise
        function at that point. */

    public static double value(double x, double y, double z) {
        double xi = Math.floor(x), yi = Math.floor(y), zi = Math.floor(z);
        double u1 = x - xi, v1 = y - yi, w1 = z - zi;
        double u2 = 1.0 - u1, v2 = 1.0 - v1, w2 = 1.0 - w1;
        int i = ((int) xi) & 0xFF, j = ((int) yi) & 0xFF, k = ((int) zi) & 0xFF;

        double unbounded =  GRAD_SCALE
                * (getWavelet(i, j, k, u1, v1, w1)
                        + getWavelet(i + 1, j, k, u2, v1, w1)
                        + getWavelet(i, j + 1, k, u1, v2, w1)
                        + getWavelet(i + 1, j + 1, k, u2, v2, w1)
                        + getWavelet(i, j, k + 1, u1, v1, w2)
                        + getWavelet(i + 1, j, k + 1, u2, v1, w2)
                        + getWavelet(i, j + 1, k + 1, u1, v2, w2) + getWavelet(
                        i + 1, j + 1, k + 1, u2, v2, w2));
        
        double bounded = Math.min(1.0, Math.max(-1.0, unbounded));
        
        return bounded;
    }

    /** Evaluate the wavelet at node (i, j, k) for the point (u, v, w). */

    private static final double getWavelet(int i, int j, int k,
            double u, double v, double w) {
        double u2 = u * u, v2 = v * v, w2 = w * w, drop;

        drop = 1.0 - u2 * (10.0 * u + (6.0 * u - 15.0) * u2);
        drop *= 1.0 - v2 * (10.0 * v + (6.0 * v - 15.0) * v2);
        drop *= 1.0 - w2 * (10.0 * w + (6.0 * w - 15.0) * w2);
        switch (grad[permute[permute[i] + j] + k]) {
        case 0:
            return drop * (u + v);
        case 1:
            return drop * (-u + v);
        case 2:
            return drop * (u - v);
        case 3:
            return drop * (-u - v);
        case 4:
            return drop * (u + w);
        case 5:
            return drop * (-u + w);
        case 6:
            return drop * (u - w);
        case 7:
            return drop * (-u - w);
        case 8:
            return drop * (v + w);
        case 9:
            return drop * (-v + w);
        case 10:
            return drop * (v - w);
        case 11:
            return drop * (-v - w);
        default:
            return 0.0; // Should never get here
        }
    }

    /** Given a point in 3D space, calculate the gradient of the scalar noise
        function at that point.  This is necessary when using noise for bump mapping. */

    public static void calcGradient(Vector3d gradient, double x, double y,
            double z) {
        double xi = Math.floor(x), yi = Math.floor(y), zi = Math
                .floor(z);
        double u1 = x - xi, v1 = y - yi, w1 = z - zi;
        double u2 = 1.0 - u1, v2 = 1.0 - v1, w2 = 1.0 - w1;
        int i = ((int) xi) & 0xFF, j = ((int) yi) & 0xFF, k = ((int) zi) & 0xFF;

        gradient.set(0.0, 0.0, 0.0);
        addWaveletGradient(gradient, i, j, k, u1, v1, w1);
        addWaveletGradient(gradient, i + 1, j, k, u2, v1, w1);
        addWaveletGradient(gradient, i, j + 1, k, u1, v2, w1);
        addWaveletGradient(gradient, i + 1, j + 1, k, u2, v2, w1);
        addWaveletGradient(gradient, i, j, k + 1, u1, v1, w2);
        addWaveletGradient(gradient, i + 1, j, k + 1, u2, v1, w2);
        addWaveletGradient(gradient, i, j + 1, k + 1, u1, v2, w2);
        addWaveletGradient(gradient, i + 1, j + 1, k + 1, u2, v2, w2);
        gradient.scale(2.0);
    }

    /** Add the gradient of the wavelet at node (i, j, k) for the point (u, v, w). */

    private static final void addWaveletGradient(Vector3d gradient, int i,
            int j, int k, double u, double v, double w) {
        double dropu, dropv, dropw, product, dot;
        double u2 = u * u, v2 = v * v, w2 = w * w;

        dropu = 1.0 - u2 * (10.0 * u + (6.0 * u - 15.0) * u2);
        dropv = 1.0 - v2 * (10.0 * v + (6.0 * v - 15.0) * v2);
        dropw = 1.0 - w2 * (10.0 * w + (6.0 * w - 15.0) * w2);
        product = dropu * dropv * dropw;
        switch (grad[permute[permute[i] + j] + k]) {
        case 0:
            dot = 30.0 * (u + v);
            gradient.x += -dropv * dropw * u2 * (u2 - 2.0 * u + 1.0)
                    * dot + product;
            gradient.y += -dropw * dropu * v2 * (v2 - 2.0 * v + 1.0)
                    * dot + product;
            gradient.z += -dropu * dropv * w2 * (w2 - 2.0 * w + 1.0)
                    * dot;
            return;
        case 1:
            dot = 30.0 * (-u + v);
            gradient.x += -dropv * dropw * u2 * (u2 - 2.0 * u + 1.0)
                    * dot - product;
            gradient.y += -dropw * dropu * v2 * (v2 - 2.0 * v + 1.0)
                    * dot + product;
            gradient.z += -dropu * dropv * w2 * (w2 - 2.0 * w + 1.0)
                    * dot;
            return;
        case 2:
            dot = 30.0 * (u - v);
            gradient.x += -dropv * dropw * u2 * (u2 - 2.0 * u + 1.0)
                    * dot + product;
            gradient.y += -dropw * dropu * v2 * (v2 - 2.0 * v + 1.0)
                    * dot - product;
            gradient.z += -dropu * dropv * w2 * (w2 - 2.0 * w + 1.0)
                    * dot;
            return;
        case 3:
            dot = -30.0 * (u + v);
            gradient.x += -dropv * dropw * u2 * (u2 - 2.0 * u + 1.0)
                    * dot - product;
            gradient.y += -dropw * dropu * v2 * (v2 - 2.0 * v + 1.0)
                    * dot - product;
            gradient.z += -dropu * dropv * w2 * (w2 - 2.0 * w + 1.0)
                    * dot;
            return;
        case 4:
            dot = 30.0 * (u + w);
            gradient.x += -dropv * dropw * u2 * (u2 - 2.0 * u + 1.0)
                    * dot + product;
            gradient.y += -dropw * dropu * v2 * (v2 - 2.0 * v + 1.0)
                    * dot;
            gradient.z += -dropu * dropv * w2 * (w2 - 2.0 * w + 1.0)
                    * dot + product;
            return;
        case 5:
            dot = 30.0 * (-u + w);
            gradient.x += -dropv * dropw * u2 * (u2 - 2.0 * u + 1.0)
                    * dot - product;
            gradient.y += -dropw * dropu * v2 * (v2 - 2.0 * v + 1.0)
                    * dot;
            gradient.z += -dropu * dropv * w2 * (w2 - 2.0 * w + 1.0)
                    * dot + product;
            return;
        case 6:
            dot = 30.0 * (u - w);
            gradient.x += -dropv * dropw * u2 * (u2 - 2.0 * u + 1.0)
                    * dot + product;
            gradient.y += -dropw * dropu * v2 * (v2 - 2.0 * v + 1.0)
                    * dot;
            gradient.z += -dropu * dropv * w2 * (w2 - 2.0 * w + 1.0)
                    * dot - product;
            return;
        case 7:
            dot = -30.0 * (u + w);
            gradient.x += -dropv * dropw * u2 * (u2 - 2.0 * u + 1.0)
                    * dot - product;
            gradient.y += -dropw * dropu * v2 * (v2 - 2.0 * v + 1.0)
                    * dot;
            gradient.z += -dropu * dropv * w2 * (w2 - 2.0 * w + 1.0)
                    * dot - product;
            return;
        case 8:
            dot = 30.0 * (v + w);
            gradient.x += -dropv * dropw * u2 * (u2 - 2.0 * u + 1.0)
                    * dot;
            gradient.y += -dropw * dropu * v2 * (v2 - 2.0 * v + 1.0)
                    * dot + product;
            gradient.z += -dropu * dropv * w2 * (w2 - 2.0 * w + 1.0)
                    * dot + product;
            return;
        case 9:
            dot = 30.0 * (-v + w);
            gradient.x += -dropv * dropw * u2 * (u2 - 2.0 * u + 1.0)
                    * dot;
            gradient.y += -dropw * dropu * v2 * (v2 - 2.0 * v + 1.0)
                    * dot - product;
            gradient.z += -dropu * dropv * w2 * (w2 - 2.0 * w + 1.0)
                    * dot + product;
            return;
        case 10:
            dot = 30.0 * (v - w);
            gradient.x += -dropv * dropw * u2 * (u2 - 2.0 * u + 1.0)
                    * dot;
            gradient.y += -dropw * dropu * v2 * (v2 - 2.0 * v + 1.0)
                    * dot + product;
            gradient.z += -dropu * dropv * w2 * (w2 - 2.0 * w + 1.0)
                    * dot - product;
            return;
        case 11:
            dot = -30.0 * (v + w);
            gradient.x += -dropv * dropw * u2 * (u2 - 2.0 * u + 1.0)
                    * dot;
            gradient.y += -dropw * dropu * v2 * (v2 - 2.0 * v + 1.0)
                    * dot - product;
            gradient.z += -dropu * dropv * w2 * (w2 - 2.0 * w + 1.0)
                    * dot - product;
            return;
        default:
            return; // Should never get here
        }
    }

    /** Given a point (x,y,z) in 3D space, set v to the value of the vector
        noise function at that point. */

    public static void calcVector(Vector3d v, double x, double y, double z) {
        double xi = Math.floor(x), yi = Math.floor(y), zi = Math.floor(z);
        double u1 = x - xi, v1 = y - yi, w1 = z - zi;
        double u2 = 1.0 - u1, v2 = 1.0 - v1, w2 = 1.0 - w1;
        int i = ((int) xi) & 0xFF, j = ((int) yi) & 0xFF, k = ((int) zi) & 0xFF;

        v.set(0.0, 0.0, 0.0);
        addVectorWavelet(v, i, j, k, u1, v1, w1);
        addVectorWavelet(v, i + 1, j, k, u2, v1, w1);
        addVectorWavelet(v, i, j + 1, k, u1, v2, w1);
        addVectorWavelet(v, i + 1, j + 1, k, u2, v2, w1);
        addVectorWavelet(v, i, j, k + 1, u1, v1, w2);
        addVectorWavelet(v, i + 1, j, k + 1, u2, v1, w2);
        addVectorWavelet(v, i, j + 1, k + 1, u1, v2, w2);
        addVectorWavelet(v, i + 1, j + 1, k + 1, u2, v2, w2);
    }

    /** Evaluate the vector wavelet at node (i, j, k) for the point (u, v, w). */

    private static final void addVectorWavelet(Vector3d vec, int i, int j,
            int k, double u, double v, double w) {
        Vector3d g = vect[permute[permute[permute[i] + j] + k]];
        double u2 = u * u, v2 = v * v, w2 = w * w, drop;

        drop = 1.0 - u2 * (10.0 * u + (6.0 * u - 15.0) * u2);
        drop *= 1.0 - v2 * (10.0 * v + (6.0 * v - 15.0) * v2);
        drop *= 1.0 - w2 * (10.0 * w + (6.0 * w - 15.0) * w2);
        vec.x += drop * g.x;
        vec.y += drop * g.y;
        vec.z += drop * g.z;
    }
}
